!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!! 
!! This module defines an atmospheric reference state.
!!
!! \n
!!
!! Define the background reference state used to subtract the leading
!! order, hydrostatically balanced terms from the flow equations. The
!! reference state is defined by the conservative variables density,
!! energy and momentum. Many additional variables are precomputed,
!! which are often used in the computations. The reference state is
!! represented by the nodal values both in the internal (element) and
!! the side quadrature points.
!!
!! \note We assume here that the reference state has zero velocity, is
!! horizontally uniform and in hydrostatic balance. This implies that,
!! once one thermodynamic variable is prescribed, all the remaining
!! variables can be derived. In particular, given the surface pressure
!! \f$p_s=p(z=0)\f$ and the temperature profile \f$T=T(z)\f$, we have
!! \f{equation}{
!!  p(z) = p_se^{\displaystyle
!!               -\frac{g}{R}\int_0^z\frac{d\zeta}{T(\zeta)}}.
!! \f}
!<----------------------------------------------------------------------
module mod_atm_refstate

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_octave_io, only: &
   mod_octave_io_initialized, &
   write_octave

 use mod_mpi_utils, only: &
   mod_mpi_utils_initialized, &
   wp_mpi, mpi_alltoallv

 use mod_base, only: &
   mod_base_initialized, &
   t_base

 use mod_grid, only: &
   mod_grid_initialized, &
   t_grid, t_ddc_grid, &
   affmap

 use mod_dgcomp_testcases, only: &
   mod_dgcomp_testcases_initialized, &
   t_phc, phc

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_atm_refstate_constructor, &
   mod_atm_refstate_destructor,  &
   mod_atm_refstate_initialized, &
   t_atm_refstate, atm_ref, &
   t_atm_refstate_e, atm_ref_e, &
   t_atm_refstate_s, atm_ref_s, &
   write_octave

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members

 !> Reference atmosphere
 type t_atm_refstate
  real(wp) :: rho   !< density
  real(wp) :: e     !< energy density \f$ J/m^3 \f$
  real(wp) :: p     !< pressure
  real(wp) :: t     !< temperature
  real(wp) :: pi    !< Exner pressure
  real(wp) :: theta !< potential temperature
  real(wp) :: phi   !< geopotential \f$gz\f$, \f$ m^2/s^2 \f$
  real(wp) :: h     !< enthalpy density \f$ J/m^3 \f$
  real(wp) :: a     !< sound speed
 end type t_atm_refstate

 !> Pointwise reference atmosphere (element points)
 !!
 !! Add information concerning the vertical derivatives.
 type, extends(t_atm_refstate) :: t_atm_refstate_e
  real(wp) :: dz_rho
  real(wp) :: dz_e
  real(wp) :: dz_p
  real(wp) :: dz_t
  real(wp) :: dz_pi
  real(wp) :: dz_theta
 end type t_atm_refstate_e

 !> Pointwise reference atmosphere (side points)
 type, extends(t_atm_refstate) :: t_atm_refstate_s
  ! nothig to add, for the time being...
 end type t_atm_refstate_s

 ! private members

! Module variables

 ! public members

 !> The reference state: first index \f$\mapsto\f$ local node, second
 !! index \f$\mapsto\f$ element
 !<
 type(t_atm_refstate)  , allocatable, protected :: atm_ref(:,:)
 type(t_atm_refstate_e), allocatable, protected :: atm_ref_e(:,:)
 type(t_atm_refstate_s), allocatable, protected :: atm_ref_s(:,:)

 logical, protected ::               &
   mod_atm_refstate_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_atm_refstate'

 interface write_octave
   module procedure write_refstate_struct
 end interface write_octave

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 !> Module constructor.
 !!
 !! This constructor allocates \c atm_ref, \c atm_ref_e, \c atm_ref_s
 !! and initializes them. Since different reference profiles require
 !! different input values, this constructor uses pointer input
 !! arguments. Only the values which are really used need to be
 !! associated.
 !!
 !! \note Currently only constant stability states are implemented,
 !! for which the temperature profile can be integrated in closed
 !! form.
 !!
 !! Using the fact that the Brunt-V&auml;is&auml;l&auml; frequency 
 !! \f{equation}{
 !!  N^2 = \frac{g}{\theta}\frac{d\theta}{dz}
 !! \f}
 !! is constant we have
 !! \f{equation}{
 !! \theta = \theta_s e^{\displaystyle \frac{N^2}{g}z},
 !! \f}
 !! and defining
 !! \f{equation}{
 !! \Delta = 1 - \frac{g^2}{c_pN^2T_s}
 !! \f}
 !! we obtain the vertical profiles
 !! \f{equation}{
 !! \pi = \pi_s\left[ e^{\displaystyle -\frac{N^2}{g}z}
 !!   + \Delta\left(1 -e^{\displaystyle -\frac{N^2}{g}z}\right)\right]
 !! \f}
 !! and
 !! \f{equation}{
 !! T = T_s\left[ 1 + \Delta\left(e^{\displaystyle \frac{N^2}{g}z}-1
 !! \right)\right].
 !! \f}
 !! Typically we have \f$\pi_s=1\f$ and \f$\theta_s=T_s\f$. A special
 !! case is when \f$\Delta=0\f$, yielding an isothermal atmosphere.
 !! Such a case is characterized by
 !! \f{equation}{
 !! T_s = \frac{g^2}{c_pN^2}.
 !! \f}
 !! and implies stable stratification.
 !!
 !! The implementation works as follows: first the
 !! Brunt-V&auml;is&auml;l&auml; frequency scale
 !! \f$f_N=\frac{N^2}{g}=\frac{1}{\theta}\frac{d\theta}{dz}\f$ and
 !! \f$T_s\f$ are determined, and then all the thermodynamical
 !! variables are obtained. The reason for using \f$f_N\f$, which is
 !! the inverse of the vertical characteristic length, is that it is
 !! well defined also for \f$g=0\f$, in which case the vertical scale
 !! is \f$+\infty\f$ and the reference profile is independent from
 !! \f$z\f$.
 !!
 !! \note To evaluate the reference atmosphere at the quadrature
 !! nodes, both internal on on the boundaries and including the
 !! derivatives of the thermodynamical variables, each variable is
 !! first projected on the FE element space by L2 projection, and then
 !! evaluated together with its derivative. Using the orthogonality of
 !! the basis we have, for a generic variable
 !! \f$u\f$,
 !! \f{equation}{
 !!  u(x)=\sum_iu_i\varphi_i(x), \qquad
 !!  u_i=\sum_lw_l\,u(x_l)\varphi_i(x_l).
 !! \f}
 !! \bug Especially for low order FE, the projection of the reference
 !! state should be consistent with the discretization used for the
 !! flow equations, in order to ensure that zero velocity
 !! configurations are preserved.
 !<
 subroutine mod_atm_refstate_constructor(grid,base,profile_type, &
              t_ref, comm,ddc_grid)
  !> grid
  type(t_grid), intent(in) :: grid
  !> FE basis
  type(t_base), intent(in) :: base
  !> Type of reference profile
  character(len=*), intent(in) :: profile_type
  !> Temperature of the isothermal reference state
  real(wp), pointer, intent(in) :: t_ref
  !> Optional arguments for ddc computations
  integer, intent(in), optional :: comm
  type(t_ddc_grid), intent(in), optional :: ddc_grid

  logical :: lnone
  integer :: ie, is, l, i, isl
  integer, allocatable :: nel4side(:)
  real(wp) :: f_n, t_s, delta, xg(grid%m,base%m), &
    z_eq(base%m), dphi_dz(base%pk,base%m), wbase(base%pk,base%m), &
    pb(base%pk,base%ms,grid%d+1), xne
  real(wp) :: theta(base%m),  pi(base%m), t(base%m), p(base%m), &
     rho(base%m), phi(base%m), e(base%m), h(base%m), a(base%m)
  character(len=1000) :: message
  ! ddc specific variables
  integer :: nsvar, id, ins, idx(base%ms), ierr
  integer, allocatable :: offse1(:), offse2(:)
  real(wp), allocatable :: senbuf(:,:,:), recbuf(:,:,:)
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.)  .or. &
       (mod_kinds_initialized.eqv..false.)     .or. &
       (mod_octave_io_initialized.eqv..false.) .or. &
       (mod_mpi_utils_initialized.eqv..false.) .or. &
       (mod_base_initialized.eqv..false.)      .or. &
       (mod_grid_initialized.eqv..false.)      .or. &
       (mod_dgcomp_testcases_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_atm_refstate_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   lnone = .false.
   profile_type_select: select case(trim(profile_type))
    case("isothermal") ! implies constant stability
     f_n = phc%gravity/(phc%cp*t_ref) ! f_N = N^2/g
     t_s = t_ref
    case("none") ! notice that few ad hoc changes will be also needed
     f_n = 0.0_wp
     t_s = 0.0_wp
     lnone = .true.
    case default
     write(message,'(A,A,A)') &
       "Unknown reference profile '",trim(profile_type),"'."
     call error(this_sub_name,this_mod_name,message)
   end select profile_type_select

   allocate( atm_ref(base%pk,grid%ne) ,                     &
     atm_ref_e(base%m,grid%ne) , atm_ref_s(base%ms,grid%ns) )

   if(f_n.gt.0.0_wp) then
     delta = 1.0_wp - phc%gravity/(phc%cp*f_n*t_s)
   else
     delta = 1.0_wp ! this will have no effect
   endif

   allocate( nel4side(grid%ns) )
   nel4side = 0
   atm_ref_s%rho   = 0.0_wp
   atm_ref_s%e     = 0.0_wp
   atm_ref_s%p     = 0.0_wp
   atm_ref_s%t     = 0.0_wp
   atm_ref_s%pi    = 0.0_wp
   atm_ref_s%theta = 0.0_wp
   atm_ref_s%phi   = 0.0_wp
   atm_ref_s%h     = 0.0_wp
   atm_ref_s%a     = 0.0_wp
   !$omp parallel do schedule(static) &
   !$omp   private( ie , xg , z_eq , theta , pi , t , p , rho , phi , &
   !$omp        e , h , a , l , dphi_dz , i , wbase , isl , pb , is ) &
   !$omp   shared( grid , base , f_n , t_s , delta , profile_type , &
   !$omp         phc , atm_ref , atm_ref_e , atm_ref_s , nel4side , &
   !$omp         lnone ) &
   !$omp   default(none)
   elem_loop_do: do ie=1,grid%ne

     ! nodes in physical space
     xg = affmap(grid%e(ie),base%xig)

     ! compute the dimensionless coordinate N^2/g * z
     z_eq = f_n * xg(grid%d,:)

     !------------------------------------------------------------------
     !1) Evaluate all the variables in the quad. nodes

     ! set the first two variables using the profiles with N=const
     theta = t_s*exp(z_eq)
     pi    = exp(-z_eq) + delta*(1.0_wp-exp(-z_eq))
     ! the next variables are set using the profiles of theta and pi
     t   = theta * pi
     p   = phc%p_s * pi**(1.0_wp/phc%kappa)
     if(lnone) then
       p   = 0.0_wp
       rho = 0.0_wp
       phi = 0.0_wp
       ! Since f_n and t_s are set to zero, the followings are implicit
       ! theta = 0
       !     t = 0
       !     e = 0
       !     h = 0
       !     a = 0
       ! On the other hand, we have to set
       pi = 0.0_wp
     else
       rho = p / (phc%rgas*t)
       phi = phc%gravity * xg(grid%d,:)
     endif
     e   = rho * ( phc%cv*t + phi )
     h   = e + p
     a   = sqrt(phc%gamma*phc%rgas*t)
     !------------------------------------------------------------------

     !------------------------------------------------------------------
     !2) For each variable, project it and use the projection to
     !compute function values and derivatives

     ! precompute vertical derivatives in physical space
     !   \nabla\phi = B^{-T}\nabla\hat{\phi}
     ! we obtain the vertical component as
     !   \partial_z phi = (B^{-1})_{:d} \cdot \nabla\hat{\phi}
     do l=1,base%m ! quadrature node
       dphi_dz(:,l) = &
         matmul( grid%e(ie)%bi(:,grid%d) , base%gradp(:,:,l) )
     enddo
     ! precompute basis times quad. weights
     do i=1,base%pk ! quadrature node
       wbase(i,:) = base%wg * base%p(i,:)
     enddo
     ! reorder the boundary values according to side ordering
     do isl=1,grid%d+1
       pb(:,:,isl) = base%pb( : , base%stab(grid%e(ie)%pi(isl),:) ,isl)
     enddo
     
     ! projection
     atm_ref(:,ie)%rho   = matmul( wbase , rho   )
     atm_ref(:,ie)%e     = matmul( wbase , e     )
     atm_ref(:,ie)%p     = matmul( wbase , p     )
     atm_ref(:,ie)%t     = matmul( wbase , t     )
     atm_ref(:,ie)%pi    = matmul( wbase , pi    )
     atm_ref(:,ie)%theta = matmul( wbase , theta )
     atm_ref(:,ie)%phi   = matmul( wbase , phi   )
     atm_ref(:,ie)%h     = matmul( wbase , h     )
     atm_ref(:,ie)%a     = matmul( wbase , a     )

     ! evaluation: elements
     atm_ref_e(:,ie)%rho   = matmul( atm_ref(:,ie)%rho   , base%p )
     atm_ref_e(:,ie)%e     = matmul( atm_ref(:,ie)%e     , base%p )
     atm_ref_e(:,ie)%p     = matmul( atm_ref(:,ie)%p     , base%p )
     atm_ref_e(:,ie)%t     = matmul( atm_ref(:,ie)%t     , base%p )
     atm_ref_e(:,ie)%pi    = matmul( atm_ref(:,ie)%pi    , base%p )
     atm_ref_e(:,ie)%theta = matmul( atm_ref(:,ie)%theta , base%p )
     atm_ref_e(:,ie)%phi   = matmul( atm_ref(:,ie)%phi   , base%p )
     atm_ref_e(:,ie)%h     = matmul( atm_ref(:,ie)%h     , base%p )
     atm_ref_e(:,ie)%a     = matmul( atm_ref(:,ie)%a     , base%p )

     atm_ref_e(:,ie)%dz_rho   = matmul( atm_ref(:,ie)%rho   , dphi_dz )
     atm_ref_e(:,ie)%dz_e     = matmul( atm_ref(:,ie)%e     , dphi_dz )
     atm_ref_e(:,ie)%dz_p     = matmul( atm_ref(:,ie)%p     , dphi_dz )
     atm_ref_e(:,ie)%dz_t     = matmul( atm_ref(:,ie)%t     , dphi_dz )
     atm_ref_e(:,ie)%dz_pi    = matmul( atm_ref(:,ie)%pi    , dphi_dz )
     atm_ref_e(:,ie)%dz_theta = matmul( atm_ref(:,ie)%theta , dphi_dz )

     ! evaluation: sides (arithmetic average of the two elements)
     do isl=1,grid%d+1
       is = grid%e(ie)%is(isl)
       !$omp critical
       atm_ref_s(:,is)%rho   = atm_ref_s(:,is)%rho   &
                        + matmul( atm_ref(:,ie)%rho   , pb(:,:,isl) )
       !$omp end critical
       !$omp critical
       atm_ref_s(:,is)%e     = atm_ref_s(:,is)%e     &
                        + matmul( atm_ref(:,ie)%e     , pb(:,:,isl) )
       !$omp end critical
       !$omp critical
       atm_ref_s(:,is)%p     = atm_ref_s(:,is)%p     &
                        + matmul( atm_ref(:,ie)%p     , pb(:,:,isl) )
       !$omp end critical
       !$omp critical
       atm_ref_s(:,is)%t     = atm_ref_s(:,is)%t     &
                        + matmul( atm_ref(:,ie)%t     , pb(:,:,isl) )
       !$omp end critical
       !$omp critical
       atm_ref_s(:,is)%pi    = atm_ref_s(:,is)%pi    &
                        + matmul( atm_ref(:,ie)%pi    , pb(:,:,isl) )
       !$omp end critical
       !$omp critical
       atm_ref_s(:,is)%theta = atm_ref_s(:,is)%theta &
                        + matmul( atm_ref(:,ie)%theta , pb(:,:,isl) )
       !$omp end critical
       !$omp critical
       atm_ref_s(:,is)%phi   = atm_ref_s(:,is)%phi   &
                        + matmul( atm_ref(:,ie)%phi   , pb(:,:,isl) )
       !$omp end critical
       !$omp critical
       atm_ref_s(:,is)%h     = atm_ref_s(:,is)%h     &
                        + matmul( atm_ref(:,ie)%h     , pb(:,:,isl) )
       !$omp end critical
       !$omp critical
       atm_ref_s(:,is)%a     = atm_ref_s(:,is)%a     &
                        + matmul( atm_ref(:,ie)%a     , pb(:,:,isl) )
       !$omp end critical
       !$omp critical
       nel4side(is) = nel4side(is) + 1
       !$omp end critical
     enddo
     !------------------------------------------------------------------

   enddo elem_loop_do
   !$omp end parallel do

   ! We need to ensure consistent values on the shared sides
   ddc_if: if(present(comm)) then
     ! 1) alloctate the buffers
     nsvar = 9 + 1 ! number of side variables (plus tag)
     allocate( senbuf(nsvar,base%ms,sum(ddc_grid%nns_id)) , &
               recbuf(nsvar,base%ms,sum(ddc_grid%nns_id)) )
     ! 2) fill the send buffer
     allocate(offse1(0:ddc_grid%nd-1)); offse1 = 0 ! global offset
     allocate(offse2(0:ddc_grid%nd-1)); offse2 = 0 ! local offset
     do i=0,ddc_grid%nd-1
       offse1(i) = sum(ddc_grid%nns_id(0:i-1))
     enddo
     fill_do: do ins=1,ddc_grid%nns
       is  = ddc_grid%ns(ins)%i
       id  = ddc_grid%ns(ins)%id
       offse2(id) = offse2(id) + 1
       l   = offse1(id) + offse2(id)
       idx = base%stab(ddc_grid%ns(ins)%p_s2s,:)
       senbuf(1,:,l) = atm_ref_s(idx,is)%rho
       senbuf(2,:,l) = atm_ref_s(idx,is)%e
       senbuf(3,:,l) = atm_ref_s(idx,is)%p
       senbuf(4,:,l) = atm_ref_s(idx,is)%t
       senbuf(5,:,l) = atm_ref_s(idx,is)%pi
       senbuf(6,:,l) = atm_ref_s(idx,is)%theta
       senbuf(7,:,l) = atm_ref_s(idx,is)%phi
       senbuf(8,:,l) = atm_ref_s(idx,is)%h
       senbuf(9,:,l) = atm_ref_s(idx,is)%a
       senbuf(nsvar,1,l) = real(ddc_grid%ns(ins)%in,wp) ! tag
     enddo fill_do
     ! 3) communication
     call mpi_alltoallv(                                           &
       senbuf, nsvar*base%ms*offse2, nsvar*base%ms*offse1, wp_mpi, &
       recbuf, nsvar*base%ms*offse2, nsvar*base%ms*offse1, wp_mpi, &
                        comm,ierr)
     ! 4) read the received data
     read_do: do ins=1,size(recbuf,3)
       is = int(recbuf(nsvar,1,ins))
       atm_ref_s(:,is)%rho   = atm_ref_s(:,is)%rho   + recbuf(1,:,ins)
       atm_ref_s(:,is)%e     = atm_ref_s(:,is)%e     + recbuf(2,:,ins)
       atm_ref_s(:,is)%p     = atm_ref_s(:,is)%p     + recbuf(3,:,ins)
       atm_ref_s(:,is)%t     = atm_ref_s(:,is)%t     + recbuf(4,:,ins)
       atm_ref_s(:,is)%pi    = atm_ref_s(:,is)%pi    + recbuf(5,:,ins)
       atm_ref_s(:,is)%theta = atm_ref_s(:,is)%theta + recbuf(6,:,ins)
       atm_ref_s(:,is)%phi   = atm_ref_s(:,is)%phi   + recbuf(7,:,ins)
       atm_ref_s(:,is)%h     = atm_ref_s(:,is)%h     + recbuf(8,:,ins)
       atm_ref_s(:,is)%a     = atm_ref_s(:,is)%a     + recbuf(9,:,ins)
       nel4side(is) = nel4side(is) + 1
     enddo read_do
     deallocate(senbuf,recbuf,offse1,offse2)
   endif ddc_if

   ! normalize side values
   !$omp parallel do schedule(static) &
   !$omp   private( is , xne ) &
   !$omp   shared( grid , nel4side , atm_ref_s ) &
   !$omp   default(none)
   side_loop_do: do is=1,grid%ns
     xne = real(nel4side(is),wp)
     atm_ref_s(:,is)%rho   = atm_ref_s(:,is)%rho   / xne
     atm_ref_s(:,is)%e     = atm_ref_s(:,is)%e     / xne
     atm_ref_s(:,is)%p     = atm_ref_s(:,is)%p     / xne
     atm_ref_s(:,is)%t     = atm_ref_s(:,is)%t     / xne
     atm_ref_s(:,is)%pi    = atm_ref_s(:,is)%pi    / xne
     atm_ref_s(:,is)%theta = atm_ref_s(:,is)%theta / xne
     atm_ref_s(:,is)%phi   = atm_ref_s(:,is)%phi   / xne
     atm_ref_s(:,is)%h     = atm_ref_s(:,is)%h     / xne
     atm_ref_s(:,is)%a     = atm_ref_s(:,is)%a     / xne
   enddo side_loop_do
   !$omp end parallel do
   deallocate( nel4side )

   mod_atm_refstate_initialized = .true.

 end subroutine mod_atm_refstate_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_atm_refstate_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_atm_refstate_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   deallocate( atm_ref , atm_ref_e , atm_ref_s )

   mod_atm_refstate_initialized = .false.
 end subroutine mod_atm_refstate_destructor

!-----------------------------------------------------------------------
 
 subroutine write_refstate_struct(rs,var_name,fu)
 ! octave output for a vector reference state
 ! Notice that the array of structures is rearranged as a structure of
 ! arrays. The output file also includes some phisical constants.
  integer, intent(in) :: fu
  type(t_atm_refstate), intent(in) :: rs(:,:)
  character(len=*), intent(in) :: var_name
 
  character(len=*), parameter :: &
    this_sub_name = 'write_refstate_struct'
 
   ! Some physical constants
   write(fu,'(a,a)')    '# name: ','phc'
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 11' ! number of fields
   ! field 01 : eradius
   write(fu,'(a)')      '# name: eradius'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(phc%eradius,'<cell-element>',fu)
   ! field 02 : omega
   write(fu,'(a)')      '# name: omega'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(phc%omega,'<cell-element>',fu)
   ! field 03 : gravity
   write(fu,'(a)')      '# name: gravity'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(phc%gravity,'<cell-element>',fu)
   ! field 04 : cp
   write(fu,'(a)')      '# name: cp'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(phc%cp,'<cell-element>',fu)
   ! field 05 : rgas
   write(fu,'(a)')      '# name: rgas'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(phc%rgas,'<cell-element>',fu)
   ! field 06 : cv
   write(fu,'(a)')      '# name: cv'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(phc%cv,'<cell-element>',fu)
   ! field 07 : gamma
   write(fu,'(a)')      '# name: gamma'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(phc%gamma,'<cell-element>',fu)
   ! field 08 : kappa
   write(fu,'(a)')      '# name: kappa'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(phc%kappa,'<cell-element>',fu)
   ! field 09 : lgas
   write(fu,'(a)')      '# name: lgas'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(phc%lgas,'<cell-element>',fu)
   ! field 10 : p_s
   write(fu,'(a)')      '# name: p_s'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(phc%p_s,'<cell-element>',fu)
   ! field 11 : p_sf
   write(fu,'(a)')      '# name: p_sf'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(phc%p_sf,'<cell-element>',fu)


   ! Now the reference state
   write(fu,'(a,a)')    '# name: ',var_name
   write(fu,'(a)')      '# type: struct'
   write(fu,'(a)')      '# length: 9' ! number of fields

   ! field 01 : rho
   write(fu,'(a)')      '# name: rho'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(rs%rho,'<cell-element>',fu)

   ! field 02 : e
   write(fu,'(a)')      '# name: e'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(rs%e,'<cell-element>',fu)

   ! field 03 : p
   write(fu,'(a)')      '# name: p'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(rs%p,'<cell-element>',fu)

   ! field 04 : t
   write(fu,'(a)')      '# name: t'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(rs%t,'<cell-element>',fu)

   ! field 05 : pi
   write(fu,'(a)')      '# name: pi'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(rs%pi,'<cell-element>',fu)

   ! field 06 : theta
   write(fu,'(a)')      '# name: theta'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(rs%theta,'<cell-element>',fu)

   ! field 07 : phi
   write(fu,'(a)')      '# name: phi'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(rs%phi,'<cell-element>',fu)

   ! field 08 : h
   write(fu,'(a)')      '# name: h'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(rs%h,'<cell-element>',fu)

   ! field 09 : a
   write(fu,'(a)')      '# name: a'
   write(fu,'(a)')      '# type: cell'
   write(fu,'(a)')      '# rows: 1'
   write(fu,'(a)')      '# columns: 1'
   call write_octave(rs%a,'<cell-element>',fu)

 end subroutine write_refstate_struct
 
!-----------------------------------------------------------------------

end module mod_atm_refstate

